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ABSTRACT 


This research designed and simulated an adaptive attitude control system for 
the Crew Equipment/Retriever (CER) during autonomous attitude hold and large 
angle or slewing maneuvers. The CER is a proposed space robot that deploys 
from the Space Station and retrieves any lost equipment or incapacitated 
astronauts. The moment of inertia tensor for the CER and acquired target is not 
known a priori. In this research, the moment of inertia tensor is estimated by a 
Kalman filter and used to update the derived linear quadratic regulator (LQR) 
and quaternion feedback regulator (QF R) control laws. Computer simulation 
results show that during attitude hold the adaptive LQR design stabilizes the CER 
and provides a more fuel efficient controller effort: as compared with a 
previously designed nonadaptive minimum time controller and a nonadaptive 
LQR design. Computer simulation results of slewing maneuvers show that the 
adaptive QFR design provides a more fuel efficient controller: as compared with 


a nonadaptive QFR design. 
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I. INTRODUCTION 


A. CERS CONCEPT DESCRIPTION 
1. CERS Origin and Purpose 

The National Aeronautics and Space Administration (NASA), Johnson 
Space Center Space Station Projects Office sent out a Request for Proposal in 
May 1987 as part of Space Station Work Package 2. This Request for Proposal, 
including an added Amendment 7, defined a requirement to provide for the 
capability to rescue an incapacitated external-vehicular activity crewman and to 
retrieve equipment inadvertently detached from the Space Station. [Ref.1: p. L- 
2-14a] 

McDonnell Douglas Astronautics Company (MDAC) responded to this 
Request for Proposal in September 1987 with a practical, low cost retriever 
concept. This concept was referred to as the Crew/Equipment Retrieval System 
(CERS). [Ref. 2: p. 1] This overall system consisted of a crew and equipment 
retriever vehicle (CER) and other Space Station based support systems. 

The overall mission of the CER is to: deploy from the Space Station, 
acquire and capture the designated target, and return to the Space Station. A 
summary of the CERS capabilities, as defined by MDAC, is listed as follows: 
[Ref. 2: p. 9] 

1. Retrieve an 850 pound target (includes 10% safety margin); 


2. Total Deployment time of 120 minutes; 


3. Retriever activated and deployed without assistance from 
an external-vehicular activity crewman; 


4. Retriever senses own attitude, range, and range to target; 


. Retriever can be remotely operated from the Space Station; 
. Retriever accommodates a worst-case separation speed of 3.5 ft/sec; 


Retriever senses and controls its own attitude with and without a target; 


. Retriever has attitude hold and three-axis translation capability. 
2. CERS Baseline Configuration 
Figure 1 shows the CER and its Space Station support systems. More 
detailed descriptions of proposed hardware and software are contained in Ref. 2 
(pp. 20-68). 
A simple representation of the CER for attitude dynamics analysis was 
developed in Ref. 3 and is shown in Figure 2. The characteristics of the baseline 


configuration developed by MDAC are listed as follows: [Ref. 2: p. 24] 


1. 850 pounds total weight; 
. Three-axis (six degrees of freedom) stabilized; 


. Remote tele-operated free flyer; 


-_ WO WN 


. Use of 24 cold Nitrogen (N29) jet thrusters rated at 1.0 Ibf each; 
5. Attitude control is accomplished by firing thrusters in pairs; 

6. Maximum control torques: 

Roll Axis 3 ft-lbf; 

Pitch Axis 3 ft-lbf; 

Yaw Axis 4 ft-lbf. 





Figure 1. CERS Major Components [Ref. 2: p. 15] 
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Figure 2. CER Baseline Configuration [Ref. 3: p. 8] 


An important feature of the CER is its ability to control its attitude. After 
acquiring a target the CER must stabilize its attitude and perform an attitude 
reorientation prior to returning to the Space Station. 

The CER's attitude control problem is very different from the norm. 
Most spacecraft are designed such that any changes in their moments of inertia 
are minimized. The control devices are, moreover, placed so as to act along the 
principal axes of the body. This positioning minimizes any gyroscopic coupling 
between torque applied to any one axis and rotation about another axis. Neither 
of these two conditions hold true for the CER after it acquires a target since the 


target can be as massive as the CER by itself. 


B. THESIS OBJECTIVES 

The overall objective of this thesis is to design adaptive attitude control laws 
for the CER with and without a target during large angle or slewing maneuvers 
and during autonomous attitude hold for all mission phases. Previous thesis 
research [Ref. 3] modeled the CER and designed time-optimal and weighted 
time-fuel optimal single-axis control systems and tested these control systems by 
applying them to several worst case target scenarios. This research analyzes the 
complete nonlinear three-axis control problem for small angle motion or attitude 
hold and large angle or slewing motion. Attitude hold pertains to attitude 
control in the presence of small disturbances while slewing motion pertains to 
attitude reorientation. The control law designs are adaptive in that a key system 
parameter, the moment of inertia of the CER and acquired target, is not known a 
priori and must be estimated. A subsidiary goal of the slewing motion control 
law is to accomplish this reorientation in an optimal fashion: an eigenaxis 


rotation. Both control laws must be able to deal with a non-diagonal moment of 


inertia tensor since after target acquisition, the thrusters no longer act along the 


principal axes of the body. 


C. THESIS ORGANIZATION 

In Chapter II, general spacecraft attitude kinematics and dynamics are 
developed. These equations of motion are then applied to the CER. The moment 
of inertia tensors for the CER with and without a target are calculated. 

Chapter III derives two control law designs. One design is developed for 
attitude hold while another design is developed for large angle or slewing 
maneuvers. Central to each control law design is the knowledge of the CER 
moment of inertia tensor. 

- In Chapter IV, an estimation scheme is developed that provides the above 
mentioned control laws with an estimate of the CER moment of inertia tensor. 
This estimation scheme is based on a rather unusual application of the Kalman 
Filter. 

Chapter V presents the computer simulation results of applying each control 
law and estimation scheme. Control system design issues and practical 
implementation details are also discussed. | 

Conclusions based on the computer simulation results are presented in 


Chapter VI. In addition, recommendations for future research are discussed. 


Il. ATTITUDE KINEMATICS AND DYNAMICS 


The equations of motion for any rotating rigid body can be divided into two 
sets: the Kinematic Equations of Motion and the Dynamic Equations of Motion. 
Kinematics studies motion without considering the forces that cause that motion. 
The Kinematic Equations of Motion are a set of first order differential equations 
that specify the time evolution of the chosen attitude parameters. The Dynamic 
Equations of Motion, meanwhile, take into account the forces that cause 
rotational motion and express the time evolution of the angular velocity of the 


rigid body. [Ref. 4: p. 510] 


A. ROTATIONAL KINEMATICS 
1. Direction Cosine Matrix 
Any general vector r can be written in terms of its magnitude and 
direction. The direction can be represented as a unit vector referenced to some 


previously defined coordinate or reference frame as shown in Figure 3. This 


unit vector is made up of components known as direction cosines: [Ref. 5: p. 9] 


1 =1f =1[(cosa)n, + (cosB)n, + (cosy)n,] (2.1) 


where r is a Scalar magnitude, and f is a unit vector whose components are 
referenced to the three orthogonal axes of the reference frame n. Note the 
following symbolic conventions used throughout this thesis: 

1. The underline bar denotes a vector, r 


2. The hat symbol denotes a unit vector, fr; 


3. The double underline bar denotes a matrix, W. 


I, =rcosa 
r, =rcosB 


ly eos) 


IT = r[(cos)n, + (cosB)n, + (cosy) ] 





Figure 3. Direction Cosines [Ref. 5: p. 9] 


Now consider a unit vector b with components along the orthogonal 
body axes of a rigid body, i.e., a spacecraft as shown in Figure 4. Let n be a 
unit vector with components along three orthogonal directions that are fixed in 
space. These two unit vectors are related by the following transformation: [Ref. 


ap. 9] 
b=Th - (2.2) 


The transformation is defined by T, the 3x3 direction cosines matrix (DCM). 


This DCM is critical to the field of spacecraft dynamics and control and Ref. 5 
addresses several of its important properties. The most important of these 
properties are summarized below: 

1. A DCM exists for any pair of orthogonal sets of three axes; 

2. The DCM is an orthogonal matrix and its inverse equals its transpose; 

3. A DCM can be built up from successive rotations about the axes. 
The last property is best explained by an example. Define a sequence of 


reference frames related by the following transformations: 


eile (2.3) 
b=T,a; (2.4) 
é= Lb. (2.5) 


and therefore the reference frames ¢ and ni are related by: 


SPACECRAFT 


n1 


ORBIT 
"FIXED" 
REFERENCE 
FRAME 


SPACECRAFT 
ORBIT 





Figure 4. Spacecraft Body Axes 
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sR Ti: (2.6) 


In other words, one overall DCM can be formed as a product of DCMs: 


tod Le (2.7) 


2. Euler Angles 

The relative orientation of two orthogonal reference frames can be 
defined in terms of three angles [Ref. 5: p. 16]. This idea was first introduced 
by Euler in the early eighteenth century and is synonymous with the idea of 
parameterizing the previously discussed nine element DCM with only three 
independent parameters. 

The classical Euler angles, for which there are twelve distinct cases, 
define an arbitrary orientation by using the successive rotational transformation 
property of the DCM. That is, a sequence of three elementary rigid right handed 
rotations about instantaneously fixed axes is used to build up an overall DCM that 
represents the transformation from one orientation to a different orientation. 
amet. 5: .p.17] 

This thesis employs the 3-2-] or yaw-pitch-roll Euler angle sequence 
since it is most commonly used in aircraft and spacecraft applications. This 
sequence is produced by initially lining up both the body axes and the fixed or 
inertial axes. A rotation about the z or number three axis is then performed, 
producing a new y and x axis . A second rotation about the new y or number 
two axis is then enacted. Finally, a rotation about the new x axis finishes the 


rotation sequence. Figure 5 depicts the formulation of this sequence. 
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Body Axes and Inertial Axes 
are initially lined up 


A rotation is then performed 
about the z or yaw axis 









ae 













Next, a rotation about the 
y or pitch axis is performed 


Finally, a rotation about the 
x or roll axis is performed 









cosy siny 0O cos8 Q -—sin®@ 1 0 0 
T,=|-siny cosy 0 |) Oa | 0 T,=|0 cos@ sing 
0 One a sin@ QO cos@ 0 -sing cos@ 








Figure 5. The 3-2-1 Euler Angles 
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Reference 5 illustrates the general process for obtaining the kinematic 
differential equations for a chosen set of Euler angles. Wertz [Ref. 4: p. 765] 
lists the kinematic equations of motion for the twelve possible Euler angle 
representations. The kinematic differential equations for the chosen 3-2-] Euler 


angle set are as follows: 


Y =(@, sin(@) + @, cos(@)) / cos(8); (2.8) 
6 =, cos(~) — @, sin(@); (2.9) 
0 = @, + (My sin(~) + W, cos(@)) tan(6); (2.10) 


where: 
1. y is the yaw angle; 
2. 9 is the pitch angle; 
3. @ is the roll angle; 


4. The vector @ is the angular velocity vector and is composed of components 
along each of the body axes. 


3. Euler's Principal Rotation Theorem 
Junkins [Ref. 5: p. 26] states that Euler is generally credited with being 


responsible for the Principal Rotation Theorem: 


A rigid body can be brought from an arbitrary initial 
orientation to an arbitrary final orientation by a single rotation of 
the body through a principle angle about a principal line; the 
principal line being a judicious axis fixed in the body and fixed in 
space. 
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This concept, displayed in Figure 6, allows the DCM to be 
parameterized in terms of the principal angle and principal line 2. In 
mathematical terms, the principal line corresponds to the eigenvector of the 
DCM: for the eigenvalue + 1. Therefore, given any DCM one can solve for the 
principal line and angle and reduce the general angular displacement to a single 
rotation about a fixed line. [Ref. 5: p. 27] 

4. Euler Parameters 

In conjunction with the Principal Rotation Theorem, Euler defined four 

parameters in terms of the principal line and principal angle. These Euler 


Parameters are as follows: 


By = cos(o/2) (2.11) 
B, = 4, sin( / 2) 
B, = £, sin(6 / 2) 

B, = 2, sin(o / 2) 


where £,, £,, £, are the components of the unit vector along the principal line 2 
and 9 is the principal angle. 

The DCM is therefore parameterized in terms of the above Euler 
Parameters and this allows the relative orientation of two orthogonal reference 
frames to be represented by four parameters. One of these parameters is 
redundant since the DCM was previously shown to be parameterized by three 
Euler angles. This redundancy manifests itself in the following constraint 


equation: 


14 





Figure 6. Euler's Principal Rotation Theorem [Ref. 5: p. 25] 
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By” as B,’ + B,” + B, =1. (2212) 


Various algorithms have been developed that determine the Euler 
parameters from a given DCM and convert back and forth between Euler angles 
and parameters. These algorithms have been used extensively in the computer 
simulation programs for this thesis and are included in the Appendix. [Ref. 5: 
pp. 31-35] 

By differentiating the inverse relationships between the Euler 
parameters and the elements of the DCM, and making a few substitutions, the 
kinematic differential equations in terms of the Euler parameters can be 


formulated as: [Ref. 5: p. 35] 


SSIs = (2:13) 


Note that throughout the literature on attitude dynamics and control, the 
Euler parameters are sometimes referred to as quaternions and are formulated as 


follows: [Ref. 4: p. 414] 


q, = 4, sin(g / 2) (2.14) 
q. = 4, sin(9 /2) 

q, = £,sin(o /2) 

q, =cos(o / 2) 
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The constraint equation (2.12) is also applicable and a conversion from one 
representation of Euler parameters to this quaternion representation allows the 


following identification: 


qs =B, (2.15) 
q; =B;, Tale 205", 


As seen above in equation (2.15), the difference between these two 
representations is rather insignificant. Note that this author has chosen to use the 
Euler parameter representation of equation (2.11) but throughout this thesis the 
terms Euler parameters and quaternions are used interchangeably. 
5. Parameterization Discussion 

The previous sections briefly demonstrate that there exist several 
choices when representing the orientation of a rigid body. Euler angles are 
easier to visualize and are more popular but suffer from one large draw back: 
the presence of mathematical singularities at certain angles. Equation (2.8), for 
example, experiences a singularity when the cosine of the pitch angle goes to 
zero. This can be a real problem in terms of numerical computations during 
computer simulation or software running control system algorithms. The use of 
Euler parameters eliminates the use of trigonometric functions and their 
singularities but they are harder to visualize. 

Therefore, the intended application should dictate the choice of 
kinematic differential equations. In this thesis, the 3-2-7 Euler angle set is used 


to define the kinematics for small angle motion during attitude hold. For large 


17 


angle or slewing motion and control, the less widely used but much more 


practical Euler parameters are used to represent the kinematics. 


B. ROTATIONAL DYNAMICS 
1. General Rigid Body 
The Eulerian Rotational Equations of Motion for a rigid body subject to 
applied control torques are well known and are typically represented as: [Ref. 5: 


pp. 49-52] 


“Glo +I "u (2.16) 


where: 
1. @ is the angular acceleration vector; 


2. I is the moment of inertia tensor for the rigid body; 


Se 


@ is the angular velocity vector; 


4. u is the vector of applied control torques; 


5. @ is the skew symmetric angular velocity matrix defined as: 
0 -@, O, 
@=/0, O —, |. (28 
®, O, 0 


2. The CER 
The dynamics of the CER can be represented by equation (2.16) 
provided that the rigid body assumption is used. The only system parameter in 
equation (2.16) that is specific to the CER is its moment of inertia tensor. This 


moment of inertia tensor will change whenever the CER captures a target. 
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Hansen [Ref. 3] calculated the moment of inertia for the CER without 
target by assuming an 850 pound total system weight symmetrically distributed 
about the center of gravity. The resulting matrix is diagonal in the CER's body 


coordinate system and has the units of slug-foot?: 


396 0 0 
lag =| 0 55 O|. (2.18) 
Diets 5 


The total moment of inertia of the CER with a target, as previously 
discussed, is not known a priori since it depends on the specific mass and 
moment of inertia of the acquired target. Hansen [Ref. 3: p. 17] utilized the 
parallel axis theorem to calculate the total moment of inertia for the CER and an 


acquired target 


Ww 


M,M 
Lora = Teer + Tarcer "Mah F: m R, (2219) 
—— rm 1 2 


where M, is the mass of the CER, M, is the mass of the target, and R, is the 
skew symmetric matrix derived from the vector r, between the center of gravity 


of the CER and the target: 


> agit amma os 
eir 0 =k (2.20) 
ait 0 


The fact that vectors can only be added together if they are expressed in the same 


coordinate system also applies to the moments of inertia in equation (2.19). Here 


A, 


again, the DCM T plays the role of transforming a moment of inertia from one 
reference frame to another reference frame as follows: 


ae = Tliarcer ig (2.21) 


where the superscript T refers to taking the transpose of the given matrix, the 
superscripts b and n refer to the reference frame, and I is the moment of inertia 
tensor. 

Hansen [Ref. 3: pp. 17-21] further utilized equations (2.18-2.21) to 
calculate five worst case target scenarios, each with an associated total moment of 
inertia tensor. This thesis used the Case Two scenario most frequently for 
control system design analysis and computer simulation. This scenario 


corresponds to an astronaut with manned maneuvering unit in the target net and 


has a total moment of inertia ( slug-foot?) of: 


112.9 24 -111.9 
loom =| 24 534.9 64 |. (2.22) 
-1119 6.4 497.6 





C. STATE VARIABLE REPRESENTATION 

The total set of equations representing the CER attitude kinematics and 
dynamics is composed of equation (2.16) and either equations (2.8-2.10) for the 
3-2-1] Euler angle representation or equation (2.13) for the Euler parameter 
representation. These equations are nonlinear in nature but for initial control 
system design and as an approximation when small angle motion is considered, 


equation (2.16) and equations (2.8-2.10) can be approximated as: 
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X= AX+Bu (2.23) 


where X = [y 8 © O, , w,] u is the applied control torques and: 
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In order to complete the state space equations, the plant defined by equations 


(2.23-2.25) must be accompanied by an output equation: 


Y=CX+Du. (2.26) 


This thesis assumes that all the states are measurable and there is no direct 
coupling between control input and the state outputs. Therefore, the following 


definitions are made: 
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Note that the C matrix is simply the identity matrix and the D matrix is simply a 


matrix of zeros. Therefore, equation (2.26) can be reduced to: 


I< 
i 
IPS 


(2.255 


Realistically, equation (2.27) should contain some added errors or noise due to 
the fact that sensors have limited accuracy and do introduce noise into any actual 
system. A more detailed discussion on sensors is contained in Chapter III and a 
related discussion on computer simulation implementation details is found in 


Chapter V. 
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II. ATTITUDE CONTROL LAW DESIGN 


The development of the attitude equations of motion in the previous chapter 
and the design of control laws or algorithms in this chapter are based on one 
fundamental assumption: attitude motion of a spacecraft can be approximately 
decoupled from its orbital motion. For the purpose of attitude control design, 
therefore, the spacecraft is almost universally considered to have only rotational 
degrees of freedom about its center of mass which is fixed to a reference frame 
moving on the orbital path. In reality, attitude and orbital dynamics are coupled 
and environmental torques produced by gravity gradient, aerodynamic, and solar 


radiation pressure depend on the spacecraft's orientation. [Ref. 6: p. 7] 


A. DESIGN PHILOSOPHY 

The goal of control system design is to cause the output variable of some 
dynamic system or process to follow a desired reference variable accurately in 
spite of changes in this reference variable, the external disturbances applied, and 
any changes in the dynamics of the process itself. Prior to beginning any control 
System design, a mathematical model of the system to be controlled is 
constructed. [Ref. 7: p. 17] The dynamic process to be controlled in this thesis 
is the attitude of the CER plus any acquired target and its equations of motion 
have been developed in the previous chapter. 

The process of regulation defines a situation in which the output variable of 
some dynamic process must follow a constant, usually zero, reference variable. 
[Ref. 7: p. 107]. The attitude control of the CER can be defined as a regulation 


process. While maneuvering or after acquiring a target, the CER's attitude must 
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remain constant despite the presence of external disturbance forces. After target 
acquisition, the CER must perform a reorientation before initiating a 
translational maneuver that will return it to the Space Station. In this 
reorientation maneuver, the reference variable is the desired orientation and the 
output variable is the current orientation. A regulation process is needed to 
reduce the orientation error, which represents the difference between the current 
and desired orientation, to zero in a timely fashion. 

This thesis designs closed-loop or feedback control systems because of their 
inherent ability to reject disturbances and errors in the model of the dynamic 
process to be controlled. This decision, by definition, requires the introduction 
of an output sensor which can introduce noise into the control system [Ref. 7: 
pp. 107-113]. The availability of quality sensors to measure angular position and 
angular velocity is assumed in this research. Wertz [Ref. 4: pp. 155-201] 
describes in detail the various types of hardware available to determine a 
spacecraft's angular position and angular velocity. The most common 
instruments used are the rate gyro and rate-integrating gyro. Since both 
instruments have a long history of operation, and are relatively inexpensive, the 


assumption that quality sensors exist is, therefore, very reasonable. 


B. SPECIFICATIONS AND CONTROL DEVICES 

The typical feedback control system is designed such that the overall dynamic 
system response meets a set of predetermined specifications. These 
specifications, more often than not, must be translated into terms more readily 
understood by the control engineer. 

One set of specifications is known as time domain specifications and include 


such information as settling time, maximum overshoot, and damping ratio. In 
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terms of the CER's dynamics, this refers to how long it takes for a reorientation 
maneuver and how much orientation angle overshoot is allowed before settling 
down to the final desired value. Typical spacecraft attitude reorientation control 
requires that no overshoot occurs and this corresponds to a damping ratio of 1.0. 
Reference 1 does not specifically address the issue of settling time for CER 
attitude reorientation maneuvers or attitude regulation in the presence of 
disturbances. However, MDAC [Ref. 2: p. 9] defines a total deployment of time 
of 120 minutes. This is certainly an upper limit for any reorientation 
maneuvers. A more realistic figure is obtained from MDAC's [Ref. 2: p. 17] 
definition of major mission phases as a function of time. In this mission phase 
sequence, ten minutes is allocated to target capture. In this thesis, 70 seconds was 
selected as a reasonable settling time to accomplish any reorientation maneuver. 
Any attitude regulation, after the CER is subject to a disturbance, must be 
accomplished in a fraction of this time. A more detailed discussion of 
specifications and performance will be given in the following sections on attitude 
hold and slewing maneuvers. 

In order to achieve control over the CER some type of control device must 
be selected. Wertz [Ref. 4: p. 201-210] discusses in detail the types of devices 
available for spacecraft control. The most widely used type of device are gas jets 
or thrusters and these are the control devices that MDAC selected in their 
baseline design of the CER. 

The choice of the control actuators is, furthermore, tied to the desired 
specifications and to what is practically available. MDAC, for example, has 
chosen 1 Ibf cold nitrogen thrusters. Wertz [Ref. 4: p. 206] mentions that gas 


jets are classified as cold gas and hot gas. Hot gas jets typically produce high 
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thrust levels (>5 N or 1.12 lbf) but rely upon a chemical reaction which must 
reach steady state. Cold gas systems produce lower thrust levels (< 1 N or 0.225 
Ibf) but do not rely upon a chemical reaction which must reach steady state. 
Cold gas systems, therefore, provide more precise control and can be used 
effectively in a pulsed mode. This thesis assumes the use of the cold gas thrusters 
selected by MDAC. These thrusters are commercially available and can operate 


in a pulsed mode. 


C. ATTITUDE HOLD 
1. Optimal Control Theory 


A linear feedback control law is defined in the following form: 


u=—-KX (3.1) 


where: 
1. u is the applied control effort; 


2. K is a gain matrix that must be determined; 


3. X is the vector containing all the state variables, assuming that they are all 
available by either measurement or estimation. 


The gain matrix K can be determined by choosing appropriate Laplace domain 
closed-loop pole locations based on some given time domain specifications. This 
method, known as pole placement, only works well for single-input-single-output 
(SISO) systems. In multi-input-multi-output (MIMO) systems, such as the CER, 


this technique does not lead to the development of a unique control law. In 
optimal control theory, the gain matrix K is determined by minimizing a 


specified performance criterion or cost function. [Ref. 8: pp. 337-338] 
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2. Linear Quadratic Regulator 

The linear quadratic regulator is an optimal control law of the form 
shown in equation (3.1). It is called linear because the control law is a linear 
function of the system states and it is called a regulator because this type of 
control law is well suited for regulation type problems. The quadratic 
description relates to the fact that the gain matrix K, of equation (3.1), is 
calculated by minimizing a quadratic integral cost function. For the continuous 
state space system described by equation (2.23), a quadratic integral cost function 


can be formulated as: 
I= J" EX" (QX(t)+ u"(Ru(t)]de (3.2) 


where Q and R are symmetric weighting matrices that must be chosen by the 
control system designer. Q penalizes deviation of the state vector X from the 
origin and R penalizes the use of too much control effort. [Ref. 8: pp. 339-340] 

During attitude hold of the CER, the goal is to reject all disturbances 
and maintain a constant attitude with zero angular velocity. If all the states are 
considered initialized at zero, then the control effort must drive all states towards 


the origin after a disturbance causes a deviation from this situation. A typical 
value for Q that will cause the position angles to go to the origin and yet limit 


the angular velocities is 
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where c is a number less than one [Ref. 8: p. 340]. 

The choice of a proper control weighting matrix R also requires 
careful consideration. Unless a cost or penalty is imposed for using too much 
control, the design that emerges may generate control signals that can not be 
achieved by the actuators or control devices. The resulting control signal then 
saturates at the maximum signal value that can be produced and this produces, in 
most cases, the fastest possible response. The fastest possible response may be 
highly desirable but the closed-loop behavior of a system in saturation may be 
quite different from the closed loop behavior predicted by a system not in 
saturation. The system may even become unstable when the control system 
saturates and because of this consequence, R should be chosen to avoid 
saturation. [Ref. 8: p. 341] 

Hansen [Ref. 3] designed a control law based on saturating the CER's 
control input. This control scheme produced a minimum time solution but it 
required the control thrusters to be either on or off; this is known in the 
literature as bang-bang control. Some of the worst case target scenarios 
produced unstable results and this caused Hansen to choose larger values and 
different locations for the CER thrusters in order to avoid an unstable situation. 


The gain matrix K, that minimizes equation (3.2), is found by solving 


the Riccati Equation. It is, in general, a time varying matrix that, given fairly 
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general conditions which apply to the CER, eventually reaches a steady state 


value. This optimum gain matrix in the steady state is given by: 


=R"B'M (3.4) 


where M satisfies the algebraic equation also known as the algebraic Riccati 


equation: 


Q=MA+A’M-MBR™B'M+Q. (3.5) 


[Ref. 8: pp. 345-346] Many software packages, including the program 
MATLAB used in this thesis, contain subroutines that calculate the steady state 
value of K for a given dynamic system and cost function. 

Attitude hold for the CER is, therefore, accomplished by using a Linear 
Quadratic Regulator to drive the states of the system to the origin after 


experiencing some external disturbances. The gain matrix K is determined by 


supplying a MATLAB subroutine with a model of the CER's dynamics and 
appropriately chosen Q and R weighting matrices. The steady state value of K 


is used for simplicity since K only varies near the final time and this situation is 


not encountered during most of the CER's mission. 


D. SLEWING MANEUVERS 

The problem of reorienting a spacecraft from one rest orientation to another 
rest orientation, although a problem of regulation, requires the formulation of a 
quite different control law. The linearized equations of motion, used for attitude 


hold design, are no longer valid. Slewing over a potentially large range of 
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orientation angles requires the use of the complete nonlinear equations of motion 
defined by equations (2.13) and (2.16). A linear control law may not be 
adequate for this task and the formulation of a nonlinear control law may be 
required. In addition, any slewing maneuver should, ideally, be an optimal 
maneuver. Optimal in this context refers to taking the shortest angular path. 
1. Ejigenaxis Rotations 

Many spacecraft attitude control systems are currently based on a 
sequence of rotational maneuvers about each control axis. This is a natural bias 
based on the popularity of Euler angles for describing rigid body orientation. 
However, the maneuver time of such successive rotations is two to three times 
longer than that of a single maneuver about the eigenaxis. This eigenaxis is the 
principal axis developed in Euler's Principal Rotation Theorem. Euler, 
moreover, proved that the principal angle 6 of equation (2.11) is always smaller 
than the algebraic sum of three successive Euler angles and represents the 
shortest angular path between two relative orientations. Therefore, a control law 
that causes a spacecraft to reorient itself by rotating about the eigenaxis will be 
executing an optimal maneuver. [Ref. 9: pp. 375-376] 

2. Quaternion Feedback Regulator 

Reference 9 develops a nonlinear control law that takes into account the 
complete nonlinear attitude equations of motion and executes an eigenaxis 
rotational maneuver. This development begins by defining a quaternion error 
which represents the attitude error between the current orientation and the 


desired or commanded orientation: 
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(3.6) 


where terms in the matrix represent the commanded orientation expressed in 
Euler parameters. A three-dimensional error vector can be formed by 


extracting the three vector components of equation (3.6) as: 


B.. 
B. =| B,. |- (3.7) 
B,. 


The complete quaternion regulator feedback control law is: 


u=@]a-Do-Kp, (3.8) 


where the first term is a nonlinear body-rate feedback term that counteracts the 
gyroscopic coupling torque found in equation (2.16), the second term is a linear 
body-rate feedback term, the third term is a linear error-quaternion or Euler 


parameter-error feedback term, and D and K are 3x3 constant gain matrices to 


be properly determined. [Ref. 9: p. 376] 


To complete the control law in equation (3.8), the matrices D and K 


must be determined. Reference 9 considered the gain selections 


ee (3,9) 
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where k and d are positive scalars, and I is the spacecraft moment of inertia 


tensor. Global stability via Lyapunov stability analysis was proved for the 


control law of equation (3.8) provided that 
K"D>0. (3.11) 


Equation (3.11) is always guaranteed with the gains defined by equations (3.9- 
3.10). 

With global stability guaranteed for the control law in equation (3.8), 
all that remains is a proper choice of k and d. Let i be a unit vector along the 


eigenaxis. The orientation is then expressed as: 


B=sin(o/2)A. (3.12) 


Substituting equation (3.8) into equation (2.16) yields the following closed-loop 


equation: 


@ = -I' Glo +I" (Gla - Do - Kp). (3.13) 


Further substitution of equations (3.9-3.10) produces: 


& =1'(-dlw — kIB) (3.14) 


= 
— 


o2 


or by rearranging terms: 


I@ + dIl@ + kI1B =0 (3315) 


When the angular rate @ is small, and when an eigenaxis rotation is assumed, the 


angular rate can be approximated as: 
O=OA. (3.16) 
Further substitution of equations (3.16) and (3.12) into equation (3.15) yields: 


($ +o +ksin(o/2))A =0. (3.17) 


Since the moment of inertia I is, by definition, a positive definite matrix and the 


unit vector 2 is non zero, then rn #0 and equation (3.17) becomes: 
b+db+ksin(o /2)=0. (3.18) 
If the sin(d/2) is approximated by 6/2, equation (3.18) is further reduced to: 
$6 +db+ko/2=0. (3.19) 


Equation (3.19) is the well-known linear second order equation where the 


damping ratio € and the natural frequency ,, satisfy: 
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d=26an; (3.20) 


k=20%, (3.21) 


Therefore, proper selection of the damping ratio and the natural frequency 
defines the positive scalars d and k. [Ref. 9: p. 377-378] 
Since this thesis has assumed a damping ratio of one, a required 


maneuver settling time is converted to a required natural frequency by: 


T, =4/Cay. (3.22) 


To account for the nonlinear effect of sin(?/2) when @ is large, equation (3.22) 


was modified by Ref. 9 as: 


T, =8/Coy. (3.23) 


Therefore, a quaternion feedback regulator control law is defined by a desired 


orientation and by a desired maneuver settling time as: 


d=2w, =16/T,; (3.24) 
=20 = 928) (3.25) 
u= G1 - (16/ T, )I@ — (128/ Ts )IB, . (3.26) 
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IV. ADAPTATION LAW DESIGN 


Knowledge of the moment of inertia tensor for the CER and acquired target 
is required by both the control laws developed in Chapter III. In the case of the 
large angle control law, imprecise knowledge of the moment of inertia tensor 
results in a non-eigenaxis slewing maneuver. An argument can also be made that 
precise knowledge of the moment of inertia tensor should limit the amount of 
thruster firings required and save propellant fuel. Therefore, the ability to 
estimate the moment of inertia tensor and provide this information to the control 
algorithms is a very beneficial addition to the previously developed control 


algorithms. 


A. ADAPTIVE CONTROL THEORY 

An adaptive controller differs from a static or ordinary controller in only 
one respect. In an adaptive controller, the controller parameters are variable 
and there is a mechanism for adjusting these parameters on-line based on signals 
available from the overall system. The two main approaches for constructing 
adaptive controllers are: the model reference method and the self-tuning 
method. A schematic representation of each of these methods is presented in 
Figure 7. (Ref. 10: p. 315] 

This thesis employs the self-tuning method. The design of an adaptive 
controller by the self-tuning method involves choosing a control law based on 
variable parameters and choosing an adaptation law for adjusting those 
parameters. The controller, therefore, couples a previously designed control law 


with an on-line parameter estimator. [Ref. 10: pp. 319-323] The previous 
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Figure 7. Adaptive Control Methods [Ref. 10: pp. 315, 320] 
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chapter developed two control laws based on knowledge of the CER moment of 
inertia tensor; the next section develops an adaptation law based on the Kalman 


filter equations. 


B. KALMAN FILTER DESIGN 
1. General Kalman Filter Equations 
The Kalman filter is an observer or state estimator for a dynamic 
process given by the following discrete state space equation: 


X(k +1) = @X(k) + A, u(k) + A, w(k) (4.1) 


— 


where: 
1. X(k) is the state vector at the present time k; 
. X(k+1) is the state vector one time step in the future; 
. @ is the discrete-time version of the A matrix given in equation (2.23); 
. u(k) is the applied control; 


w(k) is an unknown random input called the plant driving noise; 
A, is the discrete version of the B matrix given in equation (2.23); 


IAW PF WN 


A, is the random input influence matrix and is often identical to A,. 


In addition, the measurements of the system are given by: 
¥(k) = €X(k) + v(k) (4.2) 


where v(k) is arandom vector known as measurement noise. [Ref. 11: p. 27] 


oy 


The Kalman filter provides an optimal estimate of the system described 
by equation (4.1) and (4.2) by minimizing the mean square error between the 


actual states and the estimated states. The Kalman filter has the following form: 


K(k + 1k +1) = X(k + Ik) + G(k +1) ¥(k+1)- ¥(k+ 1k)| (4.3) 


where: 


1. the notation (k +1Ik) refers to the discrete value at time k +1 based on data 
accumulated through time k; 


. X(k+1lk) is the estimate of the states given data through time k; 
. ¥(k+1Ik) is the estimate of the measurements given data through time k; 


. G(k +1) is known as the Kalman filter gain; 
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. X(k+1Ik +1) is estimate of the states given data through time k +1. 
A set of recursive equations, in the proper order, that solve for the Kalman filter 


gain and equation (4.3) are: 


P(k + Ik) = @P(kIk)@" + A, WA, *; (4.4) 

G(k +1) =P(k +1Ik)CT[CP(k +11k)C™ +VJ ; (4.5) 
P(k + 1k +1) =[1-G(k +1)C]P(k + 1k); (4.6) 
X(k + 11k) = OX(klk) + A,u(k); (4.7) 

Y(k + Ik) = CX(k + 1k); (4.8) 
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X(k +k +1) = X(k + 1k) +G(k+ vy Spe Kk + k)}; (4.9) 


where: 


1. P is the covariance matrix of the estimation error; 
2. W is the covariance matrix of the plant driving noise; 
3. V is the covariance matrix of the measurement noise; 
4. 1 is an appropriately dimensioned identity matrix. 
Note that final implementation of equations (4.4-4.9) requires an initial estimate 
of the states X(0I0) and an initial estimate of the covariance matrix P(Ol0) [Ref. 
11: pp. 27-29]. The initialization of equations (4.4-4.9) is discussed in Chapter 
V. 
2. Linear Model 
For attitude hold, the linear quadratic regulator requires a model of the 
CER's dynamics and this is provided by equations (2.23-2.25). The B matrix 
defined by equation (2.24) contains the elements of the CER's inverse moment of 
inertia tensor which is the system parameter that needs to be estimated. 
To apply the previously defined Kalman filter equations, a plant 


equation similar to equation (4.1) is formulated as: 
I(k +1) =I" (k) + w(k) (4.10) 
where I"'(k) is a vector that contains the six independent elements of the inverse 


moment of inertia tensor and w(k) is random plant driving noise. By comparing 


equations (4.1) and (4.10), the following correspondences are noted: 
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Q=1y (4.11) 
1 = 9; (4.12) 


A,=1. (4.13) 


This basically assumes that the inverse moment of inertia does not change, other 
than by the addition of some random noise, from one time step until the next. 
An accompanying measurement equation is formulated from a 


linearized version of equation (2.16) as 
@=I"u. (4.14) 


Equation (4.14) can be algebraically rearranged into the following form: 


u(k) uk) u(k) 0 0 0 
@k)=| O uk) O uk) uj(k) O [I"(K)+v(k) (4.15) 
0 0 u, (k) 0 u,(k) u,(k) 


where v(k) is the combined measurement noise that results when angular 
acceleration and applied control torques are measured. By comparing equations 


(4.2) and (4.15), the following correspondence is noted: 
u,(k) u,(k) u,(k) 0O 0 0 


Cik)=| 0 uk) O- uk) uj(k) O j. (4.16) 
0 O uk) O- uk) u,(k) 
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The previously defined Kalman filter can now be used to estimate the 
inverse moment of inertia tensor of the CER plus any acquired target. This on- 
line estimation assumes that: 

1. Reasonable angular acceleration measurements or estimates are available; 
2. The applied control torques are available; 

3. The covariance matrices for plant and measurement noise can be computed; 
4 


. Initial estimates for the inverse moment of inertia and covariance matrix of 
estimation error can be computed. 


All the above assumptions are discussed further in Chapter V. 
3. Nonlinear Model 
An estimate of the moment of inertia for the CER and any acquired 
target is required by the quaternion feedback regulator in order to accomplish 
any slewing maneuvers. A plant equation similar to equation (4.10) is formed 


as: 

I(k + 1) = I(k) + w(k) CoN), 
where equations (4.11-4.13) are still valid. Since large angle motion is 
inherently nonlinear, it is only prudent to form a measurement equation from the 
original nonlinear version of equation (2.16) which has been rearranged as: 

u=@1o +Io. (4.18) 


Further algebraic manipulation results in the following form: 


u(k) = C(k)I(k) (4.19) 


4] 


where: 


Q, @,-0,0, @,+@,0, -O,0, 05;-@, @,0, 
C(k)=| @,@, @,+0,0, 05-f @, @,-@0, —@,0,{ (4.20) 
-0,0, Of -05 @,-0,0, 0,@, @,+0,@, 0, 


and the time dependence of each term is assumed. Note that while the dynamics 
of the CER are nonlinear, the dynamics of the CER's moment of inertia remain 
linear. 

The previously defined Kalman filter can now be used with this model 
to estimate the moment of inertia tensor of the CER plus any acquired target. 
This on-line estimation scheme utilizes the same assumptions listed in the 


previous section. 
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V. RESULTS 


A. SIMULATION PROGRAM AND IMPLEMENTATION 

Attitude control performance for the CER and target was evaluated through 
computer simulation. The simulation programs were written in MATLAB and 
are listed in the Appendix. The continuous-time equations of motion given in 
Chapter II and the continuous-time control algorithms given in Chapter III were 
simulated using discrete-time versions of these equations. A small sampling 
period of 75 milliseconds was chosen to ensure that the discrete equations 
accurately represented the continuous-time equations. 

The selection of thruster size and location inherently limits the amount of 
available control torques. To ensure realistic simulation results, hard limits have 
been included in the simulation programs so that control signals greater than 
those given in Chapter I are not generated. 

Plant noise has been added into the simulation programs to further increase 
the realism of the computer simulation results. Typical spacecraft disturbance 
torques are due to thruster misalignment and solar pressure. Kaplan [Ref. 12: p. 
241] lists some assumed values for disturbance torques. The maximum 
magnitudes of these assumed values are: 

1. Thruster misalignment torque: 8.5 x 10-5 N-m or 6.27 x 10-> ft-lbf; 

2. Solar Pressure torque: + 1.0 x 10-4N-m or +7.4x 10-5 ft-lbf. 
AS a worst caSe guess, plant noise was modeled as a random three dimensional 
vector with a normal or Gaussian distribution and a standard deviation equal to 


the sum of these two disturbance torques, 1.4 x 10-4 ft-lbf. 


43 


The control laws developed in Chapter II assume the availability of sensors 
to measure angular position and angular velocity. The sensors have limited 
accuracy and introduce measurement noise into the control system. Wertz [Ref. 
4: pp. 199-200] states that typical rate gyros have a resolution of less than 0.01 
degree per second and typical rate-integrating gyros have a resolution of 0.003 
degree. Reference 13 indicated that accuracies of 0.01 degree and 0.0003 degree 
per second are quite reasonable. As a worst case guess, measurement noise was 
modeled as two random three dimensional vectors with a normal or Gaussian 
distribution and standard deviations equal to 0.003 degree and 0.0001 degree per 


second. 


The linear quadratic regulator used for attitude hold required a choice of Q 
and R weighting matrices. After some trial and error, weighting matrices were 
obtained that allowed the states to be driven towards the origin in a reasonable 
amount of time (ten seconds) and also avoided saturation of the applied control 


torques: 


1 D000 0960 
0 1 0 0 06 
OO" 1 Ou 
= : mie 
Q-lo 0 000 0) (Og 
0 0 0 0.006 
000000 
2x10” 0 0 
R=| 0° 200° “Wile 62) 
0 0 2x10° 


Finally, the adaptation algorithm developed in Chapter IV required some 
computations and assumptions. Reference 13 indicated that angular 
accelerometers exist for spacecraft applications but are very expensive and not as 
common as rate gyros and rate-integrating gyros. In this thesis, angular 
acceleration measurements were generated by numerically differentiating the 


given angular velocity measurements as follows: 


@(k) = (@(k + 1) - @(k))/T, (5.3) 


where T, is the sample period or interval. Alternative methods for generating 


angular acceleration measurements might involve filtering the already available 
angular velocity measurements or using an extended Kalman filter to estimate the 
moment of inertia tensor and the angular accelerations simultaneously. The 


initial covariance of estimation error matrix P(0I0) was calculated empirically by 


computing the total moment of inertia tensor for several targets in various 


locations within the CER's net and then computing the overall mean and standard 


deviation. The covariance matrix of plant driving noise W and the covariance 


matrix of measurement noise V were computed as follows: 


W=wl; (5.4) 


V=vI1; (5.5) 


where 1 is an appropriately dimensioned identity matrix, and v and w are scalar 


values equal to the variances of the assumed plant noise and measurement noise. 
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The plant is the CER and target moment of inertia and this author assumes that it 
does not change very significantly. A small number should, therefore, be chosen 
for w, with zero an ideal choice. However, simulations indicated that a zero 
plant noise covariance matrix leads to unsatisfactory results. In this thesis, the 


following small number was used: 


w= in. (5.6) 


The measurement noise takes into account the errors and noise that are 
introduced by thruster misalignment, and measuring angular acceleration and 
angular velocity. After some trial and error via computer simulations, the value 


chosen empirically for v was: 


v=1x10". (ome 


The adaptation scheme constantly updates the control algorithms as it 
estimates the moment of inertia. The Kalman filter has some inevitable 
transients prior to convergence and computer simulation results indicated that the 
Kalman filter may even provide an estimate for the moment of inertia that is not 
physically realistic. Since the moment of inertia tensor is by definition a positive 
definite matrix (the eigenvalues are always positive for a real physical body), a 
test has been incorporated within the control and adaptation schemes: estimates 
of the moment of inertia are only passed to the control law if all the eigenvalues 


are positive. 
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Hansen [Ref. 3: p. 21] computed the total moment of inertia tensor for the 
CER and several worst case target capture scenarios. These moments of inertia, 
tabulated in Table 1, are used in the following sections to evaluate the control 


laws developed in this thesis. 


TABLE 1. TARGET CAPTURE MOMENT OF INERTIA TENSORS 
Case/Description Moment of Inertia Tensor 
(slug-foot2) 


1. No Target 


UO a 24 = 111.9 


Zo 534.9 6.4 
unit in net center not aligned with -1119 6.4 497.6 


2. Man + manned maneuvering 


CER. 
3. 850 lb point mass at net edge 


along the x-axis. 
4. 850 lb point mass at net edge 


along the y-axis. 


5. 850 lb point mass at net edge 


along the z-axis. 


6. 850 lb point mass at X=Z, Y=1 


from net center. 
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B. ATTITUDE HOLD 

Evaluation of attitude hold for the CER and acquired target was 
accomplished by initializing the position angles and angular velocities at some 
reasonable values and then allowing the control system to drive these values 
towards the origin. This thesis employed the same the initial conditions as 
Hansen [Ref. 3: p. 43]: 2.0 degrees for all position angles and 0.2 degrees per 
second for all angular velocities. Throughout the attitude hold simulations, the 
adaptive control system has been initialized with the moment of inertia defined 
by Case 1 in Table 1. 

1. Comparison of Adaptive and Nonadaptive Control 

The control system developed in Chapter III was compared with a 

nonadaptive version of the same control system for a CER and target moment of 
inertia defined by Case 2 in Table 1. The nonadaptive controller was given a 
constant value for the CER and target moment of inertia: defined by Case 1 in 


Table 1. The simulation results are listed below in Table 2. 


TABLE 2. ADAPTIVE/NONADAPTIVE SIMULATION RESULTS 


I= Case 1 


A figure of merit for any spacecraft control system design is the amount of fuel 






used. In this thesis, the sum of the absolute value of all control used for each 


design was calculated. This value is listed as fuel used in Table 2, although it is 
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actually only proportional to the amount of fuel used. Note that the settling time 
in Table 2 is the maximum settling time for all three position angles. 

Simulation plots for the nonadaptive design (I= Case 1) and the adaptive 
design are displayed in Figure 8. Both designs provide for an overall stable 
closed loop system. However, the nonadaptive design is clearly much more 
oscillatory, requires a long period for all oscillations to completely dampen out, 
and requires more control effort. The adaptive design starts out with the worst 
possible guess (the CER alone) for the moment of inertia but rapidly and 
correctly estimates the actual moment of inertia of the CER and target. The final 
estimate is very close to the actual moment of inertia defined by Case 2 in Table 


i: 


ime 7 Mess —1115 
[=| 28° S40 3.3 |. (5.8) 
sits 23 498.1 
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b. Adaptive Attitude Hold 
Figure 8. Simulation Results for Attitude Hold 


50 


2. Adaptive Control with various target scenarios 
The adaptive control system developed for attitude hold in Chapter III 
was simulated for all the target-capture moment of inertia tensors defined by 
Hansen [Ref. 3]: Cases 1-6 in Table 1. The simulation results for this adaptive 
attitude control design and the nonadaptive minimum time controller designed by 


Hansen [Ref. 3] are listed below in Table 3. 


TABLE 3. ATTITUDE HOLD SIMULATION RESULTS 


Target Case Settling Time (sec) Settling Time (sec) 





The adaptive control system was stable and exhibited reasonable settling times 
for all target cases. The minimum time control system designed by Hansen [Ref. 
3] exhibited unstable results for three of the target capture cases: as depicted by 
the infinite settling times. 

Table 4 lists the final estimate of the CER and target moment of inertia 
tensor for each of the target cases. The results are close but not exact. The 
differences may be attributed to the use of noisy measurement signals and the 


short time (10 seconds) available for system measurements. As the control 
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torque and angular acceleration measurements become very small, the estimation 
scheme is essentially provided with no new information and therefore can not 


accurately estimate the moment of inertia tensor. 


TABLE 4. MOMENT OF INERTIA ESTIMATES 





Actual Moment of Estimate of Moment 
Inertia Tensor of Inertia Tensor 
' 39.6 0 O 39.6 -0.3 0.3 
0 55 O -—0.3 55.2 0.2 
0 0 55 | 0.3 0.2 54.4 





2.4 
534.9 
6.4 






=) 
6.4 
497.6 


112 fees Sel], 5 
28 340 3.3 
ellis 33 A961 
69.3 —7.0 -179.4 
—7.0 1155.8 —-4.4 
-179.4 -4.4 1131.6 


112.9 
2.4 
Soles 




























} 989 -110.5 -110.5] | 98.7 -110.6 -109.9 
~110.5 496.1 -29.7 ~110.6 284.5 182.0 
-~110.5 -29.7 496.1 109.9 182.0 782g 


173.7 90.2 =—286.1 
—90.2 35155) 526 
—286.1 -52.6 746.8 


| yo | ere 411.4 5.0 —417.0 
0 796.4 0 50 udm —5.9 
3684 0 ~ Neen 417.9 5.9 “S209 





C. SLEWING MANEUVERS 


Evaluation of slewing maneuvers for the CER and acquired target was 


accomplished by: initializing the position angles and angular velocities at zero, 


a2 


selecting a final orientation in terms of the chosen 3-2-/] Euler Angles, and 
selecting a desired maneuver settling time. The adaptive quaternion feedback 
regulator developed in Chapters III and IV was compared with two nonadaptive 
quaternion feedback regulators that assumed a CER and target moment of inertia 
as defined by Case 1 in Table 1. The adaptive quaternion feedback regulator was 
initialized with this same moment of inertia. A desired orientation of 50 degrees 
for each position angle was selected for each controller design. The simulation 


results for a target defined as Case 2 are listed below in Table 5. 


TABLE 5. SLEWING MANEUVER SIMULATION RESULTS 





The adaptive design clearly meets the desired settling time and desired 
overshoot requirements. The first nonadaptive design appears to use more fuel 
than the adaptive design, but it is very oscillatory and does not meet the desired 
settling time. A second nonadaptive design was generated by adjusting the 
desired settling time and damping ratio, which changes the parameters k and d in 


the control law developed in Chapter III, until the actual settling time approached 
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the original desired settling time of 70 seconds. This second nonadaptive design 
meets the required settling time but the response is more oscillatory and uses 
more fuel: as compared with the adaptive design. Figures 9, 10, and 11 further 
clarify the differences between the adaptive and nonadaptive designs. 

The simulation plots of the quaternions or Euler parameters can be used to 
check for an eigenaxis rotation: which is defined by Ref. 9 as a straight line in 
any plot of the quaternions or Euler parameters. The first nonadaptive 
controller is clearly executing a noneigenaxis rotation. The second nonadaptive 
design controller and the adaptive controller are, almost exactly, executing an 
eigenaxis rotation. The adaptive controller is, therefore, the only design that: 

1. minimizes the desired position angle overshoot; 
2. meets the desired settling time; 

3. minimizes the fuel used; 

4, executes an eigenaxis rotation; 

5 


. has the ability to react to different and changing moments of inertia. 
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Figure 9. Simulation Results for Nonadaptive (Design 1) Slewing 
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Figure 10. Simulation Results for Adaptive Slewing 
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Figure 11. Simulation Results for Nonadaptive (Design 2) Slewing 
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VI. CONCLUSIONS 


A. ATTITUDE HOLD 

The adaptive attitude control system developed in this thesis to maintain 
attitude hold for the CER and an acquired target provided superior control 
during the capture of all previously defined worst case target scenarios. This 
adaptive linear quadratic regulator provided stable results with very reasonable 
settling times, and used a modest amount of fuel as compared with the previously 
designed nonadaptive minimum time control system [Ref. 3] and a nonadaptive 
linear quadratic regulator design. The additional computational burden of 
adaptive control is, therefore, compensated by a substantial improvement in 


control system performance. 


B. SLEWING MANEUVERS 

Fairly accurate knowledge of the moment of inertia tensor is essential for 
slewing maneuvers that are well damped, accomplished in a timely manner, use 
minimal fuel, and are executed as eigenaxis rotations. The adaptive quaternion 
feedback regulator developed in this thesis clearly provides this type of slewing 
maneuver control and outperforms nonadaptive quaternion feedback regulators. 
Again, the additional computational burden of adaptive control is more than 


compensated for by a substantial improvement in control system performance. 


C. FUTURE RESEARCH 
Both adaptive control system designs developed in this thesis are practical 
designs that will result in a more reliable and fuel efficient Crew/Equipment 


Retriever for the Space Station Freedom. These designs, however, can also be 
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applied to more general spacecraft attitude control problems. The adaptation 
scheme can easily be altered to account for spacecraft moments of inertia that are 
not only unknown but are subject to frequent and significant changes. In this 
respect, the adaptive linear quadratic regulator is an ideal candidate for large 
space structures of the future. Instead of limiting space structure design to 
shapes and load distributions that lead to quick approximations of the overall 
moment of inertia and minimize any changes, this adaptive control system would 
remove all these restrictions and become, in current science fiction terminology, 
an automatic attitude damping system. 

Future research could examine alternative adaptation schemes and employ 
filtering techniques to obtain angular acceleration measurements from angular 
velocity measurements, since the latter are more generally available. Additional 
computer simulations could attempt to model the CER and target as a flexible 
Structure and account for any target movement within the CER's capture net 


mechanism. 
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APPENDIX. MATLAB SIMULATION PROGRAMS 


Jo PROGRAM NAME: CERSimNL1 
%o PROGRAM AUTHOR: _ Nicholas F. Russo 
Jo CALLED PROGRAMS: _actualmoi, CERLQR, guessmoi 
%o 
% DESCRIPTION: 
Jo The below program will simulate the CER Control System for Small 
% Angle Motion or Attitude Hold. The Control Law is a Steady State Linear 
% Quadratic Regulator. However, it is an Adaptive Controller in that the 
% Inverse MOI of the CER + Target is unknown and potentially variable. A 
J% KALMAN Filter is used to estimate the Inverse MOI and then update the 
% feedback gains produced by the Linear Quadratic Regulator (LQR). It is 
% assumed that the angular accelerations are available as measurements. 
%o 
clear 
% Prompt user for first guess inverse MOI tensor 
guessmoi; 
%o 
% Define the initial state estimate for the Kalman Filter 
est_x=[lig(1,1);lig(1,2);lig(1,3);lig(2,2);lig(2,3);lig(3,3)]; 
Jo 
% Define the Q and r parameters for the LOR 
Q=diag([1,1,1,0,0,0}); % The Q matrix weights the system states 
Jo 
t=2.0e-5; % This is used to weight the control inputs used 
% This is the for used R=r*eye(3,3); 
%o 
% Define the actual inverse MOI tensor, as Ii. 
%o This is placed into the MATLAB workspace by a separate 
% subroutine. 
actualmoi; 
Jo 
T=75e-3; % The sampling time. 
B=[zeros(3,3);li]; % The B matrix. 
‘0 
t(1)=1; %Initial sampling index 
% Initial conditions in degrees & degrees/sec and in radians & radians/sec 
x, I)=(2;2;2;0.2;0.2;0.2}.*(pi/180) 
Ye 


oO 
y(:,1)=x(;,1); % Initial measurements 
% 
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% Initial Kalman Filter Parameters 
phik=eye(6); % Define the phi matrix for the plant 


Jo 
Ew=(le-8)*eye(6); % Plant Noise Covariance Matrix 
V=1.0e-4; 


Ev=V*eye(3); % Measurement Noise Covariance Matrix 

%o Prompt user for Covariance Estimation Error 

ques3=input(‘Specify the Covariance Estimation Error, P(O/O). < O > The 
Default. <1 > The Identity Matrix. <2 > A matrix of zeros.>>>'); 

Jo 

if ques3==0, 

oldP=[-0.0114 -1.14e-6 2.92e-5 2.166e-5 4.84e-7 2.95e-5 

-1.14e-6 0.0001 -2.56e-7 -1.9e-7 -4.25e-9 -2.587e-7 

2.92e-5 -2.56e-7 -0.0026 4.864e-6 1.087e-7 6.623e-6 

le-6 le-6 le-6 -.0019 le-6 le-6 

le-6 le-6 le-6 le-6 -4.246e-5 le-6 

le-6 le-6 le-6 le-6 le-6 -2.587e-3]; 

elseif ques3==1, 

oldP=eye(6); 

else, 

oldP=zeros(6,6); 

end 

% This is P(O/O), the covariance matrix of the estimation error 

Jo 

% Prompt user for number of simulation steps 

ques4=input(‘Enter the number of Simulation Steps. <0 > The Default, which is 
160. Otherwise, enter your own value (within reason....unless you have all 
day!!>>>'); 

if ques4==0, 

NUM=160; 

else, 

NUM=ques4; 

end 

Jo 

rand(‘normal'); 


0 
rd=(1.4e-4)*rand(3,NUM); % Random Plant Disturbance 
mn1=(1.745e-4)*rand(3,NUM); % Random Measurement Noise 
mn2=(5Se-6)*rand(3,NUM);% 
mn=[mn1;mn2]; 


%o 
Yo Now, simulate the CER + Target 
for N=1:NUM 
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Jo 

wskew=[0 -x(6,N) x(5,N);x(6,N) O -x(4,N);-x(5,N) x(4,N) 0]; 
Jo Skew symmetric matrix 

a=-1*]i*wskew*inv(Ii); 

iS=tan(x(2,N))*sin(x(1,N)); 

i6=cos(x(1,N))*tan(x(2,N)); 

jo=cos(x(1,N)); 

j6=-sin(x(1,N)); 

k5=sin(x(1,N))/cos(x(2,N)); 

k6=cos(x(1,N))/cos(x(2,N)); 

A=[0 0 0 1 15 16;0 0 0 0 j5 j6;0 00 0 k5 k6;0 0 O a(1,1) a(1,2) a(1,3);0 0 O a(2,1) 
a(2,2) a(2,3);0 O O a(3,1) a(3,2) a(3,3)]; 

% 

% Now, discretize the state space equations 

[phi,delul J=c2d(A,B,T); 


oO 
% Call the CERLQR function to determine the Steady State Feedback Gain 
J Matrix 
K=CERLQR(Q,r,lig); 
%o 
t(N+1)=N+1; 9% Increment time index 
% 
%o 
UC,N)=-K*y(,N); 
% Now, apply the limits for the control torques 
if abs(U(1,N)) > 4 
Ud,N)=4*sign(U(1,N)); 
else 
Ud,N)=U(1,N); 
end 
if abs(U(2,N)) > 3 
U(2,N)=3*sign(U(2,N)); 
else 
U(2,N)=U(2,N); 
end 
if abs(U(3,N)) > 3 
U(3,N)=3*sign(U(3,N)); 
else 
U(3,N)=U(3,N); 
end 
%o 
torque(:,N)=U(:,N)+rd(:,N); 
Jo The plant discrete state equations 
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x(:,N+1)=phi*x(:,N)+delu] *torque(:,N); 
%o Measurement Equation 
y(:,N+1)=x(:,N+1)+mn(:,N); 

Te 


(8) 
xd(:,N)=A*x(:,N)+B*U(:,N); 
yka(:,N)=xd(4:6,N); % This measurement equation takes the angular velocities 
% and obtains angular accelerations (like an accelerometer ). 
% 
yk(:,N)=(y(4:6,N+1)-y(4:6,N))./T; % This equation obtains angular accelerations 
% by numerical differentiation. 
Jo 
% Now, call the Kalman Filter in order to estimate the 
% actual Inverse MOI and update the LQR Gain Matrix. 


tl=U(1,N); 

t2=U(2,N); 

t3=U(3,N); 

Ck=[t1 t2 t3 00 0;0 t1 0 t2 13. 0;0 0 tl O t2 t3); 

% 

% Below are the Kalman Filter and gain equations 


% First, the Gain Equations 
newP=phik*oldP*phik' + Ew; 
G=newP*Ck'*inv(Ck*newP*Ck' + Ev); 
oldP=(eye(6) - G*Ck)*newP; 


est_xx=phik*est_x(:,N); 
est_y= Ck*est_xx; 
est_x(:,N+1)= est_xx + G*(yk(:,N)-est_y); 
%o 
I11=est_x(1,N); 
I12=est_x(2,N); 
ba 3=est_x(3,N); 
[22=est_x(4,N); 
I23=est_x(5,N); 
[33=est_x(6,N); 
IEST=[111 112 113;112 122 123;113 123 133); 
% 
% Check the inverse MOI guess for physical reality. Eigenvalues must be > 0. 
eigcheck(1:3,N)=[eig(IEST)]; 
% 
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% Test for physical reality ( All eigenvalues of IEST >0) 
checke(N)=1 *sign(eigcheck(1,N))+1*sign(eigcheck(2,N))+1 *sign(eigcheck(3,N)) 


’ 


Update the initial guess inverse MOI 


Otherwise, do not update the control law 


% End of Simulation Loop 


% Convert the time vector to actual time 

f=trk: 

Yo create a special time vector for the control input plot 

tU=0:N-1; 

tU=tU.*T; 

Jo Calculate the total fuel used (actually this is proportional to it) 
FUELUSED=sum(sum(abs(U))) 


End of Program 


PROGRAM NAME: guessmoi 
PROGRAM AUTHOR: Nicholas F. Russo 
CALLED PROGRAMS: none 


DESCRIPTION: : 

The below program will prompt user for first guess inverse MOI tensor. 
input('Welcome to the Non-Linear Simulation of the CER and Target. Please 
choose a first guess inverse MOI tensor for the LOR. <0 > The Default. 
MAN+MMuU in Net Center with CER Frame of Reference.< 1 > The CER alone. 
<2> Case Two from Hansen Thesis. *Any other number will use an average 
of all the inverse MOI tensors>>>'); 





it S=—2, 

lig=(0.0114 -0.0001 0.00256;-0.0001 0.0019 4.246e-5;0.00256 4.246e-5 
0.002587]; 

% 

elseil S==1e 

fig=inv({[39.6 0 0;0 55 0;0 0 55)); 

elseil S=—U), 

lig=[0.0138 0 0.003;0 0.0018 0;0.003 0 0.0026] 
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% 
else. 
Jig=inv({143.8 -33.6 -175.2;-33.6 646 -10.5;-175.2 -10.5 562]) 
end 

% 

Ig=inv(lig) 

% This is the guess MOI. 

% 

% End of Program 





% PROGRAM NAME: actualmoi 
% PROGRAM AUTHOR: Nicholas F. Russo 
% CALLED PROGRAMS: | none 


Jo DESCRIPTION: 
Jo This program will question the user for a CER + target 
J moment of inertia or its inverse. 


caseN=input('Define the actual Moment of Inertia by case number>>’); 
if caseN==1, 

MOI=[39.6 0 0;0 55 0;0 0 55]; 

Jo This is the CER alone. 

% 

elseif caseN==2, 

MOI=[112.9 2.4 -111.9;2.4 534.9 6.4;-111.9 6.4 497.6]; 

% MAN + MMU 

% 

elseif caseN==3, 

MOI=[69.3 0 -178.2;0 1153.8 0;-178.2 O 1124.1]; 

% 

elseif caseN==4, 

MOI=[98.9 -110.5 -110.5;-110.5 496.1 -29.7;-110.5 -29.7 496.1]; 
% 

elseif caseN==5, 

MOI=[369.5 0 -368.4;0 796.4 0;-368.4 0 466.4]; 

% 

elseif caseN==6, 

MOI=[172.7 -93.7 -282.3;-93.7 839.7 -39.8;-282.3 -39.8 732.9]; 
% 

else, 

end 

li=inv(MOI); 


%o 
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| % End of Program 


PROGRAM NAME: CERLQR 
PROGRAM AUTHOR: _ Nicholas F. Russo 
CALLED PROGRAMS: | dlqr 


DESCRIPTION: 
The below subroutine will compute the Steady State Linear 
Quadratic Regulator Control Gains for the CER using a Linear 
State Space Model. It will receive updates on the inverse 
Moment of Inertia (MOI) Tensor from another subroutine 
that will estimate this matrix using a Kalman Filter. 


function K= CERLQR(Q,r, lh) 
Here are the input parameters: 
Q is the weighting matrix for the states 
r is the scalar used to weight the control inputs 
Ii is the inverse MOI Tensor 


function K=CERLQR(Q,r,I1); 
%o 
R=r*eye(3,3); 
O 
% Now, define the state space equations and discretize. 
T=75e-3; % The sampling time. 
A=[(000100;000010;000001;000000;000000;000000]; 
% The above is the A matrix. 
%o 
B=[zeros(3,3);li]; % The B matrix. 
% Now, discretize the state space equations 
[phi,delul J=c2d(A,B,T); 
Jo 
% Now, call the function dlqr to calculate the feedback gains 
Jo 
[K,S]=dlqr(phi,delul,Q,R); 
% END of the function 


PROGRAM NAME: replotSA 
PROGRAM AUTHOR: __ Nicholas F. Russo 
CALLED PROGRAMS: _ subplot, plot 





Jo DESCRIPTION: 
Jo Plotting Program for Small Angle Motion 
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clg; 
subplot(221); 

plot(t,x(1,:),t,x(2,:),'--',t,x(3,:),'-.'), grid 

xlabel(‘Time - seconds’), ylabel(‘Position - radians’) 
title(Position Angle Histories’) 

% 

subplot(222); 

plot(t,x(4,:),t,x(5,:),'--',t,x(6,:), '-.'), grid 

xlabel(‘Time - seconds'),ylabel('Angular velocity - rad/sec’) 
title( Angular Velocity Histories’) 

% 

subplot(223); 
plot(tU,U(1,:),tU,U0(2,:),’--",tU,UG3,:),’-.'), grid 
xlabel(‘Time - seconds'),ylabel(‘Control Torque - foot-lbs') 
title(‘Control Torque Histories’) 

% 

subplot(224); 

J Plot fuel usage plot 
plot(tU,sum(abs(U))),grid,title(‘Fuel Usage’) 

xlabel(‘Time - seconds’), ylabel(‘Magnitude’); 

pause; 

% Prompt user to see if additional plots/information are desired. 


ques5=input(‘Would you like to see the final estimate of the Inverse Moment of 
Inertia Matrix?? If the answer is yes then type 0, otherwise type 1 .>>>'); 

if quesS==0, 

hi 


subplot(221); 
plot(t,est_x(1,:)),grid 

title(’ Inverse Moment of Inertia Estimate- I] 1') 
xlabel(‘Time - seconds’) 
subplot(222); 
plot(t,est_x(2,:)),grid,title(‘112') 
xlabel('Time - seconds’) 
subplot(223); 

plot(t,est_x(3,:)), grid,title(‘I13') 
xlabel(‘Time - seconds’) 
subplot(224); 
plot(t,est_x(4,:)), grid, title(I22') 
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xlabel('Time - seconds’),pause; 


plot(t,est_x(5,:)),grid,xlabel(Time - seconds’) 

title’ Inverse Moment of Inertia Estimate - I23') 

subplot(222); 

plot(t,est_x(6,:)),grid,xlabel(‘Time - seconds’) 

title(I33'), pause; 

cle; 

plot(eigcheck(1,:)),grid,title(First Eigenvalue of the Estimated Inverse MOI 
Tensor'),pause 

plot(eigcheck(2,:)),grid,title(Second Eigenvalue’),pause 
plot(eigcheck(3,:)),grid,title(Third Eigenvalue’),pause 


End of Program 


PROGRAM NAME: CERSlewK 1 
PROGRAM AUTHOR: Nicholas F. Russo 


CALLED PROGRAMS: _ guessmoi, actualmoi, EULERa2p, 
EULERp2a321 


DESCRIPTION: 
The below program will simulate the CER Control System for Large 
Angle (Slewing) Motion. The Control Law is a Quaternion Feedback 
Regulator. However, it is an Adaptive Controller in that the MOI 
of the CER + Target is unknown and potentially variable. A KALMAN 
Filter is used to estimate the MOI and then update the feedback gains 
produced by the Quaternion Feedback Regulator. This KALMAN Filter 
assumes that angular acceleration measurements are available. 

clear 

Prompt user for first guess inverse MOI tensor 

Jo 

guessmoi; 

% Define the initial state estimate for the Kalman Filter 

est_x=[Ig(1,1);Ig(1,2);Ig(1,3);1g(2,2);1g(2,3);1g(3,3)]; 

% 





Jo Define the actual MOT tensor. 

%o This is placed into the MATLAB workspace by a separate subroutine. 
actualmoi; 

Jo 

T=75e-3; % The sampling time. 
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t(1)=1; %Initial sampling index 
% Initial orientation 
anginit=[0 0 0]; 

J 


‘O 
xEULER(1:3,1)=[anginit(1);anginit(2);anginit(3)]; 
Jo Prompt the user for the desired orientation 
angc=input('Specify the final or desired orientation in terms of 3-2-1 Euler 
Angles, in degrees. Enter as follows, [120 120 120]>>>>>’); 
% 
% Now, convert each of these orientations to Euler Parameters 
Jo Initial Conditions 
binit=EULERa2p(2,0,anginit(1),anginit(2),anginit(3)); 
beta(:, 1 )=[binit(1);binit(2);binit(3 );binit(4)]; 
w(:,1)=[0;0;0]; 
% 
bemd= EULERa2p(2,0,angc(1),angce(2),angc(3)); 
Jo Compute the commanded quaternion matrix 
bOc=bemd(1); 
blc=bcmd(2); 
b2c=bemd(3); 


b3c=bcmd(4); 

BCMD=[b0c ble b2c b3c;-blc bOc b3c -b2c;-b2c -b3c bOc blce;-b3c b2c -blc 
bOc]; 

%o 

J Prompt the user for the values used to compute the K and D weighting 

Jo matrices used in the Quaternion Feedback Regulator. 


settle=input(‘Please enter the desired Settling Time in seconds>>>>'); 
oO 

runtime=input(‘Enter desired simulation run time (sec)>>>>'); 

actual=runtime/T; 

NUMeE=actual+0.2*actual; 

NUM=round(NUM); 

Jo 

% The two weighting matrices, K and D will now be calculated as follows: 

zeta=1.0 

Jo 

omegaN=8/(zeta* settle); 

dscalar=2*zeta*omegaN; 

kscalar=2*omegaN“2; 

Jo 

rand(‘normal'); 

%o 

rd=(1.4e-4)*rand(3,NUM); % Random Plant Disturbance 


69 


mn1=((1.745e-4)/3)*rand(4,NUM); % Random Measurement Noise 
mn2=((5e-6)/3)*rand(3,NUM);% 

Jo 

% Initial Kalman Filter Parameters 

phik=eye(6); % Define the phi matrix for the plant 

Jo Prompt user for Kalman Filter Design Plant Noise 
Wnoise=input(‘Specify the scalar multiplier for the Kalman Filter Plant Noise 
Matrix. <0> Default Value (TRUST ME). <1>No NOISE. < 2 > Choose 
your own value.>>>>>'); 

if Wnoise==0, 

W=le-8; % The scalar multiplier of the below matrix 

Jo 

elseif Wnoise==1, 

W=0; 

else, 

W=input(‘Specify your own value.>>>>>'); 

end 

Ew=W*eye(6); % Plant Noise Covariance Matrix 

V=1.0e-04; 

Ev=V*eye(3); % Measurement Noise Covariance Matrix 

Jo Prompt user again, for Covariance Estimation Error 
ques3=input(‘Specify the Covariance Estimation Error, P(0/0). < 0 > The 
Default. <1> The Identity Matrix. <2 > A matrix of zeros.>>>'); 
Jo 

if ques3==0, 

oldP=inv({-0.0114 -1.14e-6 2.92e-5 2.166e-5 4.84e-7 2.95e-5 
-1.14e-6 0.0001 -2.56e-7 -1.9e-7 -4.25e-9 -2.587e-7 

2.92e-5 -2.56e-7 -0.0026 4.864e-6 1.087e-7 6.623e-6 

le-6 le-6 le-6 -.0019 1e-6 le-6 

le-6 le-6 le-6 le-6 -4.246e-5 le-6 

le-6 le-6 le-6 le-6 le-6 -2.587e-3]); 

elseif ques3==1, 

oldP=eye(6); 

else, 

oldP=zeros(6,6); 

end 

% This is P(O/O), the covariance matrix of the estimation error 

Jo 

tclockO=clock; 

% Initialize omega with noise 

wN(:,1)=w(:,1); 

Jo Now, simulate the CER + Target 

for N=1:NUM 
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% 

wN(:,N+1)=w(:,N)+mn2(:,N); % angular velocity plus measurement noise. 
% 

betaN(:,N)=beta(:,N)+mn1(:,N); % Position measurements plus noise. 
wskewN=[0 -wN(3,N) wN(2,N);wN(3,N) 0 -wN(1,N);-wN(2,N) wN(1,N) 0]; 
% 

wskew=[0 -w(3,N) w(2,N);w(3,N) 0 -w(1,N);-w(2,N) w(1,N) 0]; 

% 

Gskew=.5*[0 -w(1,N) -w(2,N) -w(3,N);w(1,N) 0 w(3,N) -w(2,N);w(2,N) 
-w(3,N) 0 w(1,N);wG,N) w(2,N) -w(1,N) 0}; 

Ye 


(a) 

t(N+1)=N+1; 9% Increment time index 

Jo 

Jo 

%o Define the error quaternion for the feedback equation 

be=BCMD*betaN(:,N); 

qe=[be(2);be(3);be(4)]; 

Jo Define the weighting matrices K and D. Note that they are updated by the 

Jo Kalman Filter's Estimate of the MOI. 

K=kscalar*lIg; 

Dq=dscalar*Ig; 

U(:,N)=wskewN*Ig*wN(:,N) - Dq*wN(:,N) - K*qe; 

J Now, apply the limits for the control torques 

if abs(U(1,N)) > 4 

Ud,N)=4*sign(U(1,N)); 

else 

Ud,N)=U(1,N); 

end 

if abs(U(2,N)) > 3 

U(2,N)=3* sign(U(2,N)); 

else 

U2,N)=U(2,N); 

end 

if abs(U(3,N)) > 3 

U(,N)=3* sign(U(3,N)); 

else 

U3,N)=U(3,N); 

end 

Jo 

torque(:,N)=U(:,N)+rd(:,N); % Apply the random disturbance torques 
oO 

Jo Discrete Nonlinear Kinematic and Dynamic Equations 

w(:,N+1)=w(:,N)+T*(li* torque(:,N)-li*wskew*MOI*w(:,N)); 
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beta(:,N+1)=beta(:,N)+T*Gskew*beta(:,N); 
% The real angular acceleration 
wdotreal(:,N)=li*torque(:,N)-litwskew*MOI*w(:,N); 
% Angular acceleration from numerical differentiation 
wdot(:,N)=(wN(:,N+1)-wN(:,N))./T; 
Jo 
% Now, call the Kalman Filter in order to estimate the 
% actual MOI and update the Control Law. 
Io 
Jo 
Ck=[wdot(1,N) (wdot(2,N)-wN(1,N)*wN(3,N)) (WN(1,N)*wN(2,N)+wdot(3,N)) 
(-wN(2,N)*wN(3,N)) (WN(2,N)42-wN(3,N)42) (WN(2,N)*wN(3,N)) 
wN(1,N)*wN(3,N) (wdot(1,N)+wN(2,N)*wN(3,N)) (WN(3@,N)42-wN(1,N)42) 
wdot(2,N) (wdot(3,N)-wN(1,N)*wN(2,N)) -wN(1,N)*wN(3,N) 
-wN(2,N)*wN(1,N) (WN(1,N)42-wN(2,N)42) (wdot(1,N)-wN(2,N)*wN(3,N)) 
wN(1,N)*wN(2,N) (wdot(2,N)+wN(1,N)*wN(3,N)) wdot(3,N)]; 
To 
% Below are the Kalman Filter and gain equations 
Jo 
% First, the Gain Equations 
newP=phik*oldP*phik' + Ew; 
G=newP*Ck'*inv(Ck*newP*Ck' + Ev); 
oldP=(eye(6) - G*Ck)*newP; 
To 
est_xx=phik*est_x(:,N); 
est_y= Ck*est_xx; 
est_x(:,N+1)= est_xx + G*(torque(:,N)-est_y); 
%o 
I11=est_x(1,N); 
li@=est_ x3); 
113=est_x(3,N); 
I22=est_x(4,N); 
Ii23=estex(5,N); 
133=est_x(6,N); 
IEST=[I11 112 1133112 122 123;113 123 133]; 
Jo 
% Check the inverse MOI guess for physical reality. Eigenvalues must be > 0. 
eigcheck(1:3,N)=[eig(IEST)]; 
v/s 


oO 
checke(N)=1 *sign(eigcheck(1,N))+1 *sign(eigcheck(2,N))+1 *sign(eigcheck(3,N)) 


% 
if checke(N)==3, 
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% Update the initial guess MOI 


Otherwise, do not update the control law 


Convert the Euler Parameters to Euler Angles 
xEULER(:,N+1)=EULER p2a321 (beta(:,N+1)); 
end % End of Simulation Loop 
%o 
% Convert the time vector to actual time 
tt." 
% create a special time vector for the control input plot 
tU=0:N-1; 
tU=U.*T; 
etime(clock,tclock0)/60 
%o Calculate the total fuel used (actually this is proportional to it) 
FUELUSED=sum(sum(abs(U))) 


End of CER Slewing Simulation Program 


PROGRAM NAME: replot 
PROGRAM AUTHOR: _ Nicholas F. Russo 
CALLED PROGRAMS: _ plot, subplot 


DESCRIPTION: 
Plotting program for previously calculated data 


xlabel(‘Time - seconds’), ylabel(’Angle - degrees’) 
title(Z or Yaw Angle History’) 

%o 

subplot(222); 

plot(t,xEULER(2,:)), grid 

xlabel(‘Time - seconds’), ylabel(‘Angle - degrees’) 
title’ Y or Pitch Angle History’) 

%o 

subplot(223); 

plot(t,xEULER(3,:)), grid 

xlabel(‘Time - seconds’), ylabel(‘Angle - degrees’) 
title(X or Roll Angle History’) 


B 


Jo 

subplot(224); 

plot(t,w(1,:),t,w(2,:),'--',t,w(3,:),'-.'), grid 

xlabel(‘Time - seconds’), ylabel(Angular velocity - rad/sec’) 
title Angular Velocity History’) 

pause; 

clg; 

%o 

subplot(221); 

plot(beta(2,:),beta(3,:)), grid 

titleQuaternion Plot- q2 vs. ql’),xlabel(‘q1'), ylabel(‘q2’) 
%o 

subplot(222); 

plot(beta(2,:),beta(4,:)), grid 

title(‘(Quaternion Plot- q3 vs. q1'),xlabel(‘q1'),ylabel(‘q3’) 
% 

subplot(223); 

plot(beta(3,:),beta(4,:)), grid 

title(Quaternion Plot- q3 vs. q2’),xlabel(‘q2'), ylabel(‘q3’) 
%o 

subplot(224); 
plot(tU,U(1,:),tU,U(2,:),'--’,tU,0(3,:),’-.') grid 
title(‘Control Torque History’) 

xlabel(‘Time - seconds'),ylabel(‘Control Input/Torque - foot-lbs’) 
pause; 

clg; 

%o 

% Prompt user to see if additional plots/information are desired. 
ques5=input('Would you like to see the final estimate of the Inverse Moment of 
Inertia Matrix?? <O> YES. <1>NO>>>'); 

if quesS==0, 

MOI 

pause; 

MOlestimate=IEST 

pause; 

invMOI=Ili 

pause; 

invMOlestimate=inv(IEST) 

pause; 

MOleigenvalues=eig(MOI) 

pause; 

MOlesteigenvalues=eig(MOlestimate) 

pause; 
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subplot(221); 

plot(t,est_x(1,:)),grid 

title(’Moment of Inertia Estimate - I11') 

xlabel('Time - seconds’) 

subplot(222); 

plot(t,est_x(2,:)), grid, title(T12’) 

xlabel('Time - seconds’) 

subplot(223); 

plot(t,est_x(3,:)),grid,title(‘I13’) 

xlabel('Time - seconds’) 

subplot(224); 

plot(t,est_x(4,:)),grid,title(122') 

xlabel(‘Time - seconds'),pause; 

clg; 

subplot(221); 

plot(t,est_x(5,:)),grid,xlabel(‘Time - seconds’) 

title’ Moment of Inertia Estimate - 123’) 
subplot(222); 

plot(t,est_x(6,:)),grid,xlabel(‘Time - seconds’) 

title(I33') 

subplot(223); 

Jo Plot fuel usage plot 

plot(sum(abs(U))),grid,title(Fuel Usage’),pause; 
oO 


clg; 


‘O 
plot(eigcheck(1,:)),grid,title(First Eigenvalue of the Estimated MOI 
Tensor'),pause 
plot(eigcheck(2,:)),grid,title’Second Eigenvalue'),pause 
plot(eigcheck(3,:)),grid,title(‘Third Eigenvalue'),pause 

oO 

else, 

end 

Jo 

% End of Program 


PROGRAM NAME: EULERa2p 
PROGRAM AUTHOR: Nicholas F. Russo 


CALLED PROGRAMS: | none 





Jo DESCRIPTION: 
Jo The below function will compute the four Euler parameters 
% given a sequence of three Euler Angles. A 3-2-1 or a 3-1-2 
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Jo rotation sequence is allowed. 


% 

% function E=-EULERa2p(e,r,al ,a2,a3); 

% 

%o Here are the required input arguments: 

% Is the sequence 3-2-1 (e=2) or 3-1-2 (e=1)? 

% Are the three angles in radians ? (Yes: r=1, No: r=0) 
% The three Euler Angles (al,a2,a3) 

% 


% A vector is returned that contains b0,b1,b2,b3. 
function E=-EULERa2p(e,r,al ,a2,a3); 

%o 

% Convert the angles to radians if necessary 
i= 

%o 

al=al*pi/180; 

a2=a2*pi/180; 

a3=a3*pi/180; 

% 

else 

end 

Jo Define the R matrices 

RO=[1 00 0;0 100;0010;000 1); 

R1=[0 -1 00;1000;000 1;00 -1 0); 

R2=[0 0 -1 0;0 0 0 -1;1 00 0;0 1 0 0); 

R3=[0 00 -1;00 1 0;0 -1 0 0;1 00 O]; 

Jo 

% Test to see which rotation sequence is desired and proceed. 
Jo 

ife==1 % This is a 3-1-2 sequence 
R312=(cos(a3/2)*RO 
+sin(a3/2)*R2)*(cos(a2/2)*RO+sin(a2/2)*R1)*(cos(al/2)*RO+sin(al/2)*R3); 
E=R312*[1;0;0;0); 

else 

% This is a 3-2-1 sequence. 
R321=(cos(a3/2)*RO+sin(a3/2)*R1)*(cos(a2/2)*RO+sin(a2/2)*R2)*(cos(al/2)*RO 
+sin(al/2)*R3); 

E=R321*{[1;0;0;0]; 

end 

% 

% End of function 


%o PROGRAM NAME: EULERp2a321 
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% PROGRAM AUTHOR: _ Nicholas F. Russo 
% CALLED PROGRAMS: | none 


% DESCRIPTION: 
% The below function will calculate a sequence of three Euler Angles 


% from given Euler Parameters. A 3-2-1 rotation sequence is 
% allowed. 

%o 

% function E=-EULERp2a321(bvector) 

% 

% Here is the required input argument: 

% The 4 Euler Parameters: b0,b1,b2,b3 

Jo Enter these parameters as a vector. For example, 

Jo bvector=[b0;b1;b2;b3] 


% A vector is returned that contains the three angles in degrees. 

function E>-EULERp2a321(bvector); 

% 

bO=bvector(1); 

bl=bvector(2); 

b2=bvector(3); 

b3=bvector(4); 

% First, compute the Direction Cosine Matrix (DCM) from the Euler 
Parameters. 

Cp=[b042+b142-b2%2-b342 2*(b1*b2+b0*b3) 2*(b1*b3-bO*b2);2*(b1 *b2- 
bO*b3) bO*2-b142+b2%2-b342 2*(b2*b3+b0*b1);2*(b1*b3+b0*b2) 2*(b2*b3- 
bO*b1) bO42-b142-b24%2+b32]; 


%o 

Jo Calculate the middle angle, beta 
a2=atan2(-Cp(1,3),((Cp(1,1))42+(Cp(1,2))%2)*.5); 
Jo Test this angle and calculate the other angles 


a2d=a2*(180/pi); 

a2d=round(a2d); 

if a2d==90, 

al=0;% This is the angle alpha or the rotation about the z-axis 
a3=atan2(Cp(2,1),Cp(2,2)); 

v/ 


elseif a2d==-90, 

al=0; 

a3=-atan2(Cp(2,1),Cp(2,2)); 

% 

else, 
al=atan2((Cp(1,2)/cos(a2)),(Cp(1,1)/cos(a2))); 
a3=atan2((Cp(2,3)/cos(a2)),(Cp(3,3)/cos(a2))); 


vad 


end 

% 
E=[al;a2;a3]; 
E=(180/pi)*E; 
Ye 


End of function 


PROGRAM NAME: Inew 
PROGRAM AUTHOR: _ Nicholas F. Russo 
CALLED PROGRAMS: | none 


DESCRIPTION: 
The below function will calculate the new Moment of Inertia (MOT) 
Tensor for the CER plus the Target. 
The Target's MOI can be in the CER ref frame or its own; 
as long as a Direction Cosine Matrix (DCM) is available 
function I=Inew(Icer,Itar.M1,M2,r2,ref,dcm); 


Here are the required input arguments: 
MOI of CER (Icer) 
MOI of Target (Itar) 
ref (If it equals O then Itar is in the CER ref frame) 
DCM from the Target ref frame to the CER ref frame (dcm) 
Mass of the CER (M1) 
Mass of the Target (M2) 
The vector between the Center of Mass of the CER 
and the target mass in Cartesian Coord. (12) 


function I=Inew(Icer,Itar,M1,M2,r2,ref,dcm); 

%o First, test to see if the Target MOI is in the CER ref frame 
Mi tet== 

Itar = Itar; 

else 

Itar= dcm*Itar*dcm’; 

end 

% 

% Define the elements of the R2 Matrix 

R2=[0 -r2(3) r2(2);r2(3) O -r2(1);-12(2) 12(1) 0); 
MF=M1*M2/(M1+M2); % The Mass Factor multiplying R2 
% 

I= Icer + Itar - MF*R2*R2; 

% End of Function 
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% PROGRAM NAME: DCM 
Jo PROGRAM AUTHOR: Nicholas F. Russo 
Jo CALLED PROGRAMS: | none 


% DESCRIPTION: 

% The below function will calculate the Direction Cosine Matrix 

% (DCM) for a 3 Euler Angle rotation sequence. A 3-2-1 (Yaw, 
% Pitch, Roll) or a 3-1-2 (Yaw, Roll, Pitch) is allowed. 


% function D=DCMCe,r,al,a2,a3); 

Jo 

% Here are the required input arguments: 

% Is the sequence 3-2-1 (e=2) or 3-1-2 (e=1) ? 

% The three Euler Angles (al,a2,a3), in radians. 
Jo Are the Euler Angles in radians (Yes: r=1, No: r=0) 
function D=DCM¢Ce,r,al ,a2,a3); 

Jo 

% Convert the angles to radians if necessary 

if t= 

al=al*pi/180; 

a2=a2*pi/180; 

a3=a3*pi/180; 

else 

end 

% First, set up the individual rotation matrices. 
C3=[cos(al) sin(al) 0;-sin(al) cos(al) 0;0 0 1); 

% This corresponds to a rotation about the # 3 axis 

Jo 

C2a=[cos(a2) 0 -sin(a2);0 1 0;sin(a2) 0 cos(a2)] ;% #2 axis 
% 


Oo 
C2b=[cos(a3) 0 -sin(a3);0 1 0;sin(a3) 0 cos(a3)]; % # 2 axis 
Je 


oO 
Cla=[1 0 0;0 cos(a3) sin(a3);0 -sin(a3) cos(a3)]; % # 1 axis 
Jo 
C1b=[1 0 0;0 cos(a2) sin(a2);0 -sin(a2) cos(a2)];  % # 1 axis 
Jo 
Jo 
% Test to see which rotation sequence is desired and proceed. 
Jo 


if e==1 
D= C26* Cltb*G3: 
else 


= Ora C2a7C 3: 


1 


D= Cla*C2a*C3; 

end 

% 

% End of this function 
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